package math2;

/**
 * Polynomials: [a, b] -> R
 * 
 * @author hbui
 */
public class PolynomialRToR extends FunctionRnToR {

	/** Coefficients of the polynomial */
	private double[] k;

	/**
	 * Constructs the polynomial.
	 * 
	 * @param a
	 *            lower bound of domain
	 * @param b
	 *            upper bound of domain
	 * @param k
	 *            vector of coefficients of polynomials ordered in the ascending
	 *            powers (input of variable length) *
	 */
	public PolynomialRToR(double... k) {
		int n = k.length;
		int t = 0;
		// tracks down the first nonzero highest-degree coefficient
		for (int i = n - 1; i >= 0; i--) {
			if (k[i] != 0) {
				t = i;
				break;
			}
		}
		this.k = new double[t + 1];
		for (int i = 0; i <= t; i++) {
			this.k[i] = k[i];
		}
	}

	@Override
	public FunctionRnToR add(double alpha, double beta, FunctionRnToR v) {
		if (v instanceof PolynomialRToR) {
			PolynomialRToR pv = (PolynomialRToR) v;
			double[] coefs = new double[Math
					.max(this.getP() + 1, pv.getP() + 1)];
			for (int i = 0; i < coefs.length; i++) {
				coefs[i] = alpha * this.getCoefficient(i) + beta
						* pv.getCoefficient(i);
			}
			if (!allZeros(coefs)) {
				return new PolynomialRToR(coefs);
			}
			return new NullFunctionRnToR();
		}
		return super.add(alpha, beta, v);
	}

	/**
	 * Tests if the vector of coefficients is a null vector
	 * 
	 * @param coefs
	 *            vector of coefficients
	 * @return true if the vector of coefficients is a null vector and false
	 *         otherwise
	 */
	private boolean allZeros(double[] coefs) {
		for (int i = 0; i < coefs.length; i++) {
			if (coefs[i] != 0) {
				return false;
			}
		}
		return true;
	}

	@Override
	public FunctionRnToR multiply(double alpha) {
		int n = k.length;
		if (alpha != 0) {
			double[] coefs = new double[n];
			for (int i = 0; i < n; i++) {
				coefs[i] = alpha * this.k[i];
			}
			return new PolynomialRToR(coefs);
		}
		return new NullFunctionRnToR();
	}

	@Override
	public double[] gradientAt(double... x) {
		int n = this.getP();
		if (n > 0) {
			double sum = n * this.k[n];
			for (int i = n - 1; i >= 1; i--) {
				sum = x[0] * sum + i * this.k[i];
			}
			return new double[] { sum };
		}
		return new double[1];
	}

	@Override
	public int getP() {
		if (k.length == 1) {
			return 0;
		}
		return k.length - 1;
	}

	@Override
	public double valueAt(double... x) {
		// implement the Horner's scheme
		int n = this.getP();
		double sum = this.k[n];
		// if n = 0, the loop does not operate
		for (int i = n - 1; i >= 0; i--) {
			sum = x[0] * sum + this.k[i];
		}
		return sum;
	}

	public PolynomialRToR derivative() {
		int n = this.getP();
		if (n != 0) {
			double[] coefs = new double[n];
			for (int i = n - 1; i >= 0; i--) {
				coefs[i] = (i + 1) * this.getCoefficient(i + 1);
			}
			return new PolynomialRToR(coefs);
		}
		return new PolynomialRToR(0);
	}

	/**
	 * Returns the array of coefficients of the polynomial.
	 * 
	 * @return the vector of coefficients
	 */
	public double[] getCoefficients() {
		return this.k;
	}

	/**
	 * Returns the coefficients of the polynomial with index idx.
	 * 
	 * @param idx
	 *            the index of coefficient w.r.t power of variable x
	 * @return the value of the coefficient
	 */
	public double getCoefficient(int idx) {
		if (idx < this.k.length) {
			return this.k[idx];
		}
		return 0;
	}

	@Override
	public FunctionRnToR multiply(double alpha, FunctionRnToR v) {
		if (alpha == 0) {
			return new NullFunctionRnToR();
		}
		if (v instanceof PolynomialRToR) {
			PolynomialRToR w = (PolynomialRToR) v;
			int pv = this.getP();
			int pw = w.getP();
			double[] c = new double[pv + pw + 1];
			for (int j = 0; j <= pv + pw; j++) {
				for (int i = 0; i <= j; i++) {
					if (i <= pv && (j - i) <= pw) {
						c[j] += this.k[i] * w.k[j - i];
					}
				}
				c[j] = alpha * c[j];
			}
			return new PolynomialRToR(c);
		}
		return super.multiply(alpha, v);
	}

	/**
	 * divide this polynomial to (x-c) using Horner algorithm. Note that it will
	 * always return the divisor even when f(c)!=0 (i.e polynomial is not
	 * divided by (x-c))
	 * 
	 * @param c
	 * @return the divisor (polynomial)
	 */
	public PolynomialRToR divide(double c) {
		int n = this.k.length;
		if (n == 1) {
			throw new Error("cannot divide constant to a polynomial");
		} else {
			double[] coeff = new double[n - 1];
			coeff[n - 2] = k[n - 1];
			for (int i = n - 2; i > 0; i--) {
				coeff[i - 1] = c * coeff[i] + this.k[i];
			}
			return new PolynomialRToR(coeff);
		}
	}

}